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The reduced ID Poisson-Nernst-Planck (PNP) model of artificial nano- 
pores in the presence of a permanent charge on the channel wall is studied. 
More specifically, we consider the limit where the channel length exceed 
much the Debye screening length and channel's charge is sufficiently small. 
Ion transport is described by the nonequillibrium steady-state solution of 
the PNP system within a singular perturbation treatment. The quantities, 
1/A - the ratio of the Debye length to a characteristic length scale and e 
- the scaled intrinsic charge density, serve as the singular and the regular 
perturbation parameters, respectively. The role of the boundary conditions 
is discussed. A comparison between numerics and the analytical results of 
the singular perturbation theory is presented. 



PACS numbers: 05.60.Cd, 05.40. Jc, 81.07.De 
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1. Introduction 

In many physical situations an exact solution of the full problem could 
not be obtained, so that various asymptotic and perturbative techniques 
must be used. Moreover, in many cases, even regular perturbation analysis 
is useless, and singular perturbation methods must be applied. 

An example of such a case constitutes a boundary-layer problem. It is 
a ID differential-equation-boundary-value problem on the unit interval for 
which the highest derivative of the differential equation is multiplied by a 
small parameter 5: 

5y"{x) + c{x)y'{x) + d{x)y{x) = f{x), (1) 

with the boundary conditions 

y(0)=a, y{l) = b. (2) 

This boundary- value problem is singular because in the limit (5^0 one 
of the solutions abruptly disappears and the limiting solution is not able 
to satisfy the two boundary conditions in ([2]). This comes from the fact 
that by setting 6 = 0, the order of the differential equation is reduced by 
1. One of the method to solve the boundary- value problem dH) and ^ 
is to look for solutions in the form of series expansion in powers of 5. In 
order to construct an uniformly and globally valid solution, the interval 
< a; < 1 is decomposed into two kinds of regions, an outer region, in 
which the solution varies slowly as a function of x, and an inner region or 
boundary-layer region, in which the solution varies rapidly as a function of 
X. A boundary-layer region is a narrow region those thickness is typically of 
order 6 or some power of 6 [Ij. One of such boundary value problem is the 
ID steady-state PNP system with non-vanishing permanent surface charge 
on the walls of a nanopore. 

Siwy and co-workers reported [21 [3] that ion transport in nanopores with 
asymmetric fixed charge distributions is characterized by such interesting 
phenomena as ion current fluctuations, rectification, and pumping. For 
this reason, we analyse the boundary value problem of the one-dimensional 
steady-state PNP system with non-vanishing permanent surface charge. 
Since the conical geometry is experimentally relevant, we consider an uni- 
formly charged conical pore. 

The flow of ions through the nanopore caused by an externally applied 
electric field is analyzed by means of the Nernst-Planck equations together 
with the Poisson equation, in a self-consistent manner. We are interested in 
ion transport phenomena occuring in a very long channel i.e. in the limit 
where the channel length exceed much the channel radii. It justifies the 
approximation of the channel as a one-dimensional object [Ij. 



3 



The aim of the paper is to show an apphcation of the singular per- 
turbation theory to such one-dimensional PNP systems where the following 
quantities - the ratio of the Debye length to the channel length L, denoted 
as 1/A ~ and the channel surface charge a - serve as the perturbation 

parameters. The system can be viewed as a singularly perturbed problem 
in 1/A and a regularly perturbed in e that denotes a dimensionless scaled 
surface charge density. In the long channel limit, we analyzed the leading 
term in 1/A, while e is considered as the regular expansion parameter. In 
our approach, we make use of above-described method to construct an uni- 
formly and globally valid solution by calculating separately the outer and 
inner solutions and matching them then in a smooth manner. 



2. Ion Transport 

2.1. One dimensional Poisson-Nernst -Planck (PNP) 

In this paper we study ion transport through a long conical nanopore 
of the length L = 1 (scaling of the z-coordinate by the channel length L, 
whereas the radial coordinate r is scaled by the small opening radius R{0)) 
with uniformly charged wall, cf. Fig. [TJ The ion flow through the charged 
conical nanopore could be driven by the ion concentration gradients and by 
the electric field modeled together by means of the electrodiffusion equation. 
Diffusion in orthogonal direction of the pore is confined to a channel of 
variable area of t[B?{z) (see Fig. [T|). The electric field inside the pore is in 
turn governed by applied field and by the ion concentrations through the 
Poisson equation. We are interested in the non-equillibrium steady-state ion 
flux and electrical current which can persist either due to applied voltage 
or due to a concentration gradient. 

One can reduce the problem from 3D to ID assuming instantaneous 
equillibration in the transverse direction 

c(z, t) ^ c{z, r, t)A{z), ^{z, t) « $(z, r, i), (3) 

where 

A{z) = ttR\z), R{z) = (1 + 7z) (4) 

with 7 = -R(l) — 1 denoting a scaled slope of the cone's radius, with the 
scaling of the radial coordinate r by -R(O). 

The 3D PNP equations for steady state ion currents /k+ and /qj- thus 
become in the ID approximation (cf. Ref. [1] and Appendix A): 
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Fig. 1. Schema of the section of the conical channel. The channel's wall is uniformly 
charged with surface charge density a. The local radius of the cone i?(z), is given 
by Eq. H 



(5) 



d _ , , _ , A^(z) _ , , 2 dRiz) 

-^K+ = :r^K+(^) +ck+(2)— J ^K+(^)^T^— J — , 

dz dz dz 

d Miz) _ . . 2 dRjz) 

Ic\- = -r^c\-{z) - cci^{z)— (^c\~{z)——— — , 

dz dz R[z) dz 

together with the reduced Poisson equation 
d^<^{z) 2 dR{z)d^{z) . e ,2 1 / / X 



(6) 



(7) 



R{z) dz dz R{z) 7rR'^{z 

where 

1_ l eoe^kBTR^joj _ 

is the scaled Debye screening length, and the effective dimensionless surface 
charge density reads 

. = (8) 

ecoNA ^ ^ 

In Eq. dZD, c denotes dimensionless concentration (see the definition given 

by Eq. ([TOl)), Co = [eoe^fcBr/(2e2NAl03cbuik)] ^^^^^ (in meters) is the Debye 
length in bulk (see Ref. [H |5] ) , eo the dielectrical constant of vacuum, ^ 
80 the relative dielectrical constant of water, /cb the Boltzmann constant, 
T the temperature, e the elementary charge, A^a the Avogadro number, 
Cbuik the bulk concentration of ions, and cq is an arbitrary reference ID- 
concentration (in 1/6.023 x 10"^"^ moles/m). Note that the coordinate z is 
also dimensionless in the above equations (measured in units of L). 
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2.2. Rigorous Boundary Conditions 
The boundary conditions are 

Ck+(0) = 7rcK+,L, Cci-(O) = vrcci-,L, 

Ck+(1) =7r(l+7)2cK+,R, cci-(l) = ^(l + 7)'cci-,R, (9) 

$(0)=0, ^(1) = ^R, 

wherein 

q,{L,R} = 10'cbulk,{L,R}^'(0)/co, i = K+, Cr (10) 

and c^mi]^ |L Rj. denote the bulk concentrations of ions on the left and right 
sides (in moles), respectively. The electro-neutrality condition yields: Ck+,l = 
Cci- L = cl, and Ck+,r = Cci-,r = cr. Furthermore, the difference of dimen- 
sionless potentials across the nanopore is related to the applied voltage U 
(in units of Volts) by $(0) - ^>(1) = eU/{kBT), yielding U = -kBT<P{l)/e. 

3. Singular perturbation study 

3.1. Perturbation parameters 

The coresponding system contains two small parameters. The first one 
1/A is related to the ratio of the Debye length to a characteristic length scale. 
Assuming that A ^ oo we search for an aproximate solution of equations 
([5]) and ([6]) in the form 

^•(2) = $(°)(Z) + ^$W(Z) + ..., 

A 

CK+(^)=cSJi(^) + ^C«(^) + ..., (11) 

The second small parameter corresponds to the re-scaled intrinsic charge 
density e (see in Eq. [8|). Furthermore, we represent ^^^\z) , {z) , c^^^_ {z) 
as a series in e 

cl>«(z) = c^«(z)-ecl>g(z) + ..., 
c-« (z) = - (i)(^) + • • • ' (12) 
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where i = 0,1,.... Finally, a singularly perturbated boundary - value 
problem in two independent perturbation parameters 1/A and e is obtained. 
We search for the solution letting 1/A 0, while e is considered as the 
expansion parameter. Note that the perturbation problem is singular in 
1/A and regular in e. 



3.2. The outer expansion. 

These series are obtained by holding z fixed and letting 1/A ^ 0. In 
zeroth order in 1/A, the problem reduces to 

(4°l(z) - c^^liz)) /ttR^z) = (13) 

reflecting the charge neutrality in the interior, and 

r{0) _ d (0) ^R'iz)M ^ d^^^){z) JO) 



-^^i-^c-(.)-2^c-(.) + ^c-(.). 



Adding and substracting the current formulas, we obtain 



-2 



CI- dz V"K+v-/ ' -CI 



R{z) "cr^ 'J dz 



z 



(14) 



(4°|(z) + cgl(.)) - -^eTTRiz), (15) 



where -I^l+I^^^. = = const, /K++^cr = ^^^^ = The both sorts 

of ions have equal diffusion coefficients (see Ref. [0]), i.e. = -Dq- = D, 
thus J^^^ yields an approximation to the negative of the total mass current, 
whereas /^''^ approximates the total electric current (see also Ref. [7]), 

..i^(.)^^-.(0) + -^(40)(,)+,g)(,)) 



^R{z) 



(4l(.) + cgl(.)) 



(16) 



(4i(.) 



(0) , A d<^(o\z) 



cr 



dz 



+7rei?'(z). 
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Then one gets 

(4°lw+c(«u.)) (j(«)-A(o<ji(.)+cSl(.)) 



dz 

+2|M(4°i(.)+cSl(.))) = .ei?(.)(/(0)-.ei?'(^^^ ^''^ 

One cannot find exact solution to this differential equation. It is solved 
perturbatively in e with tlic solution approximated as 

(4°) (.) + cgl (.)) = (.) = 4% (.) - ecgV) (^) + . . . 

(18) 

$w(z) = $[jj(z)-6$[;5(z)+..., 

where the upper index denotes the order of expansion in 1/A, and the lower 
one in e. 

3.3. The boundary layer solutions. 

Let 5(1/ \) be the width of the boundary layer which is a function of 
1/A. The problem is rescaled near z = by setting ( = z/6{l/\). It has 
been found that S{l/X) = 1/A and for z in the boundary layer, = 0(1). 
Thus, we introduce the stretched coordinate 

C = Az, (19) 

where ( G [0, +oo) and express the concentrations and the electric potential 
in terms of this coordinate as 

CK+(C/A;l/A)=p(C; 1/A), 

cci-(C/A;l/A)=n(C;l/A), (20) 
$(C/A; 1/A) = </>(C;1/A). 

In the right boundary layer, i.e. near 2; = 1 we set the another stretched 
variable 

X = A(z-l), (21) 
where x G (— oo,0] and define analogously 

ck+(x/A + 1;1/A) =p(x;l/A), 

Cci-(X/A + I;l/A)=n(x;l/A), (22) 
$(X/A + I;l/A) = ,^(x; 1/A). 
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Each of these functions is written as an asymptotic series, namely, 
p(C; l/A) = (^) + 1^(1) (C) + . . . , 1/A) = (x) + \p^'^ (x) + • • • , 

(C; 1/A) = n(°) (C) + T"^'^ (C) + ■ ■ ■ , n(x; 1/A) = n(°) (x) + ^n^^) (x) + ■ ■ ■ , 



n 



A 



<^(C; 1/A) = (C) + (C) + • • • , kx; 1/A) = (x) + 14>^'^ (x) + . . . , 



and 

i?(C)^l, i?(x)^l + 7- 



(23) 
(24) 



The upper indixes correspond to the order of the expansion in 1/A. The 
hmit expansions is obtained by holding and x fixed and letting 1/A ^ 0. 

In turn, each of the p^^\C) , n^^\C) , (t>^'^\C) functions is written as an 
asymptotic series in e, namely, 

(C) = pgi (C) - (C) + . . . , P^'^ (X) = p1S (X) - (X) + . . . , 

n^'\0 = ngj (C) - enf^ (C) + • ■ ■ , n^°\x) = ngj (x) - ehf^ (x) + • • • , 

4>^'^ (C) = (0 - (0 + . . . , 4>^'^ (X) = 0gi (X) - (X) + . . . . 

(25) 

The lower indixes refer to the order of the series expansion in e. In the 
leading order functions in 1/A the transport equations read 

= y^^Hc) + ^p(°)(c), = ^pi^Hx) + '-^P^'Hx), 

dC dC dx dx 

= |n(0)(C) - ^n(0)(C), = ±ni^\x) - ^^(o)(x), 



(26) 

d^^(°HC) 

dC ' 

(27) 



together with the Poisson equation 

l(p^'HO-n^'\Q)]+e 



9 



and obeys the following boundary conditions 





TTCL, 




= 0, 


p[S;(0)=vr(l+7)2cR, 




= 0, 




TTCL, 




= 0, 


n[Sj(0) = ^(l + 7)2cR„ 




= 0, 


<i(o) = 


= 0, 


4>) 


= 0, 






= 0, 
(28) 



with i = 1, 2, . . .. As a next step, we integrate the left boundary layer of the 
Nernst-Planck equations (j26p 

P^'\o =pSS(o)e-^(oi(c) |i + ,^(;)(c) + (-0g)(c) + (c))} + 

(29) 

After the substitution {p, n, (p} — > {p, n, (/>} in the equations (|29|) . we obtain 
the right boundary layer solutions. Then, we put these expressions for the 
concentrations into the Poisson equation (j27p . 

3.4- Matching procedure. 

Now a matching of these two representations is needed. So we choose 
constants to make left (right) boundary layer solutions and outer expansion, 
respectively, coincide for each order of e (as 1/A ^ 0) in some intermediate 
zone between the left (right) boundary layer and the outer region, respec- 
tively. Furthemore, in the left intermediate zone an intermediate variable 
Za is introduced, where 

Zct 

and < a < 1 such as 

lim — 



lim — 

Za 

The above-described procedure is also applied to find an intermediate vari- 
able {zp — 1 = (1/A)^~^x) ill the right intermediate zone. The ensuing 



(1/A)~"z, 
(1/A)i-"C, 



oo when — ^ 0, 
A 

> when -- ^ 0. 
A 



(30) 



(31) 
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'za,Zi3 fixed and 1/A - 


-0) 










= p(0)(_ 


-oo), 


cgl(0)=nW(oo), 




= n(0)(_ 


-oo), 


$(o)(o) = 0(o)(oo), 


$(0)(1) 




-oo), 



(32) 




and 

— — (oo) = 0, — - — (-oo)=0. ['id} 
dC dx 

Adding and substracting the formulas (fT3|) and (jC.ip . (fT3|) and ()C.6p . re- 
spectively, we obtain 



4S{K+,ci-}(0) = ^4S^(0)±<! ' (34) 



where the upper sign refers to ions and the lower one to CI . To find 
the values of (p^'^^oo), we make use of 

i (^pio) (oo) - n(0) (oo)) + e = 0. (35) 

Taking into account the matching conditions (j32p together with (j34p we get 

conditions for unknown constants J^*^^ , C^^-j ,E^^^ , where i = 0, 1, 2, . . . 

The above procedure allows us to find matching conditions between the 
outer expansion and the right boundary layer as well. 



3.5. The Poisson equation in the boundary layers. 

Now we consider the Poisson equation (j27p within the boundary layer, 
and integrate them following the Ref. j^7|. We are able to find analytical 
solutions of (j27p up to the second order in e. The solutions of the Poisson 
equation (j27p depend on the boundary conditions. For ck+ l = lj 

</.[;](C) = (see also Ref. [7J), <A[;j(C) = ^ (e'^^ - l), and <Ag(C) = 0. 

Similar results hold for the solutions in the right boundary layer. If 
Ck+,r = Ccr,R, then (f)'-^(x) = ^i? (as in Ref. [7J), 

^[;](x) = ^-^(ev^>^-l),and^[°)(x)=0. 
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3.6. Uniformly valid solutions. 

To obtain an approximation that is valid uniformly on [0, 1], we add the 
boundary and outer appoximations and subtract their common limit in the 
intermediate zone. Thus, we find: 

•i>(2) = <i>j;j(z) - e{^.j;!(Az)+<i>j;jw+0j;j(A(i-2)) 



T 27rcR(l + 7)'<^S?](A(1 - z)) ^ vr(2 + 7)| + ©(e^), (37) 
where the upper sign refers to K"*" ions and the lower one to Cl~. 

4. The Donnan boundary conditions 

Approximatively, one can impose the so-called Donnan equilibrium bound- 
ary conditions at the channel ends [8j, reading: 



(38) 



^Don. N _ 3, 1 , q (^{L,R})C0 

^ ^ ^ ' 10%i?(2|L,R})^Cbulk(^{L,R}) 

where ^{l,r} = ^ x \ ^) ^'^'^ i = K+, Cr. Here, i?(2;{L^m) 

-f<(.^;{L,R}j 

denotes the nanopore radii at the left and right ends of the pore, respectively. 

Let us expand the Donnan equilibrium conditions in terms of power 
of e and afterwards, using the definition of Cj |l r}, i = K+, Cl~ (see Sec. 
Rigorous Boundary Conditions), we find our dimensionless functions at both 
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channel ends 

^2, 



vr 9 vr 4 

H e 5-^ + 

8Ci 128c3i?(z|L,R})2 



2Cii?(z|L,R}) 
1 3 3 . 



48c3i2(z|L,R})3 1280cfi?(z|L3|)5 



where the upper sign in cP°" refers to potassium ions and the lower 

one to chloride, i = K+, Cl~, = 1, and zr = (1 + 7). Note that the Don- 
nan equilibrium conditions (Eq. (20)-(21) in Ref. [5]) have been obtained 
under the assumption of local electroneutrality at the channel borders. 
Theorem 1 

Let 1/A — > 0. If C ^ +00, then the left boundary layer solutions 
p^^\C,)^n^^\C,),(j)^^\C,) given by an asymptotic series [2S\) tend to the Don- 
nan asymptotic series in z^. 

One can formulate analogue of above theorem for right boundary layer 
solutions. 

Proof: We first observe that the solution of the following transport equa- 
tion 

-/f = ^/'^(0 + ^P^'\C), Ce[o,+oo), (41) 

= ^-^°HC)-^n(0)(C), (42) 

when I^n} ~^ ^ ^^^^ equilibrium condition) tends to the Boltzmann distri- 
bution: 

p{0)(^) ^ p{o)(o)e-{0(«'(C)-</''°)(o))^ (43) 
^{0)(^) ^ ^(o)(o)e(*'»'(C)~^(")(o)). (44) 

It gives us the Donnan equilibrium condition for a given potential drop. 
On the other hand, it is known from the matching procedure that p^^\C): 
n^^\C)j 'P^^HO when ^ —>■ +00 must tend to the outer expansions c[^|(0), 
£7_(0), $(0)(0). This means that for ( +00 the condition of local elec- 
troneutrality is fulfilled. And then we end up with series given by Eqs. (j40p . 
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In turn, the outer solutions c^+{z), c'^^_{z), <I>(°)(z) (given by ([37|) and 
([36}) ) tend to the Donnan asymptotic series (00]) for z = 2:{l,r}) respectively. 
Thus, on the basis of Theorem 1 and the matching procedure, one can 
formulate 

Theorem 2 

//1/A 0, then the outer expansions c^\{z),c^^^^{z),^^^\z) approxi- 
mate the solutions of the ID-PNP system (Eq. and ^) with the Donnan 
boundary conditions satisfying the local electroneutrality at the boundaries. 

Proof: 

The boundary layer solutions calculated with the Donnan boundary con- 
ditions (j40p (after inversion to normal coordinate z) are restricted to the 
points at the boundaries approximated by these series (00]) . 

The above-presented theorem was already demonstrated in Ref. [9], in 
which a singular perturbation expansion was also used to derive the Donnan 
potential. We recall this finding, however, in a different manner. One can 
conclude that in the limit where the channel length exceeds much the Debye 
screening length, the Donnan jumps at the channel borders can be safely 
used instead of the rigorous treatment of the the boundary layers of finite 
width. 

Moreover, on the basis of the matching conditions given by (|32p and 
Theorem 1, we notice that unknown constants (present in the outer ap- 
proximation), in particular the total flux J^^^ and the electric current 
in both discussed cases of boundary conditions constitute the same series in 
e. 

5. Results 

5.1. Perturbation vs. numerics 

We shall present a comparison of the analytical results calculated by 
the above-described singular perturbation method with numerical solutions. 
The ID PNP system (equations dS])-®), with boundary conditions (equa- 
tions in ([11)), is integrated numerically by making use of a collocating method 
with adaptive meshing [10]. We choose the following set of parameters: 
R{0) = 3 nm, i?(l) = 6.616 nm, L = 200 nm, 7 = 1.2055, the room tem- 
perature (T = 298 K), the relative dielectric constant of water = 80, and 
the diffusion coefficients of the ions = D^^- = 2 x 10^ nm^/s [B]. For 
such a channel, the value of the first perturbation parameter 1/A is about 
10"^ The square of this parameter is sufficiently small. However, this is 
not enough to provide a good approximation to the PNP system with a 
varying value of the surface charge density which is not that small in the 
reality. To test the limits of regular perturbation method with respect to 
the second expansion parameter e, we consider the two different cases (a) 
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a = —0.02 e/nm^, and (b) a = —0.1 e/nm^. The first one corresponding to 
e = 0.12 is well within the perturbative treatment. For the second one with 
e = 0.6 the perturbative treatment is expected to fail, but it might work 
occasionally. 

5.1.1. Uniformly valid approximation 

The electric potential and concentration profiles for both the perturba- 
tion and the numerical solution are shown in Fig. [2l Fig. [3l and Fig. HI It is 



200 




200 



Fig. 2. (Color online) Potential profile ^{z) for cl = cr =0.1 M. Calculations are 
done for $(0) = $(L) = V (a) and $(0) V, ^ 1 V (b). Solid line and 

squares: a = —0.02 e/nm^, dashed line and circles: a = —0.1 e/nm^. Symbols 
in all cases stands for numerical solution, lines represent perturbation theory (first 
order in e). Insets depict closer look into left and right boundary layers. 
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clearly shown that the agreement with the analytical result (I36p . and (I37p 
is getting worser with the increasing value of a. Moreover, the discrepan- 
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Fig. 3. (Color online) Concentration profile C}<i+{z) for cl = cr = 0.1 M. Calcula- 
tions are done for $(0) = $(L) = V (a) and $(0) = V, $(L) = 1 V (b). Solid 
line and squares: a = —0.02 e/nm^, dashed line and circles: cr ~ —0.1 e/nm^. Sym- 
bols in all cases stands for numerical solution, lines represent perturbation theory 
(first order in e). Insets depict closer look into left and right boundary layers. 



cies between theory and numerics are more significant for growing absolute 
value of voltage. Figures [3][1] illustrate the situation in which the pertur- 
bation approximation starts to miscarry for the charge density a = —0.1 
e/nm^ in the non-equillibrium situation when the transmembrane voltage 
reaches the value of 1 V. Note, in that case the perturbation theory also 
provides unphysical negative concentrations. 

However, to assess the accuracy of power series in e, we shall systemati- 
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Fig. 4. (Color online) Concentration profile Cqi-{z) for cl = cr = 0.1 M. Calcula- 
tions are done for $(0) = V (a) and $(0) = V, $(L) = 1 V (b). Solid 
line and squares: cr = —0.02 e/nm^, dashed line and circles: a — —0.1 e/nm^. 
Symbols in all cases stands for numerical solution, lines represent perturbation 
theory (first order in e). In sets depict closer look into left and right boundary 
layers. 



cally compare the approximations including higher orders of this parameter 
with numerical solutions. Such a process for the uniformly valid solutions 
of singular perturbation theory (including boundary layer solutions) is pos- 
sible only up to second order in e. Nevertheless, on the basis of Theorem 
2, one can note that the outer expansions approximate well the solutions 
of the ID-PNP problem with the Donnan boundary conditions. One can 
obtain these perturbative expansions for any order in e. 
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5.1.2. Outer approximation 

In the non-equillibrium situation, we examine the electric potential pro- 
file for e = 0.6 (Fig. [2j)), and the concentration profile for e = 0.12 (Fig. [8)3). 
Both approximations (in the first order of e) seem to show not too large 
discrepancy from the numerical solution. However, if we plot the same ex- 
pansion in Eq. ()36p with only one more term in the perturbation series, 
we notice a considerable deviation from the numerical results one when the 
transmembrane voltage becomes sufficiently large: C/ = +lV, orC/ = — IV 
(Fig. [5|). This confirms that for the charge density a = —0.1 the perturba- 
tion method fails totally. However, for a smaller charge density a = —0.02 
the perturbation theory works as demonstrated in Fig. [Gj where the pertur- 
bation expansion was calculated up to the fourth order in e. One can notice 
that the agreement with numerics is getting better with growing order of 
the perturbative expansion. Then, it yields very accurate approximations 
in all aspects: the potential and concentration profiles, the total fiux, and 
the electric current as well. 

5.2. Rigorous b. c. vs. the Donnan b. c. 

In Fig. [71 we compare the rigorous boundary layer solution with one 
given by the Donnan potential jump approximation. The both agree well 
except from the behavior in the boundary layer, where the "Donnan" solu- 
tion makes naturally a jump because it corresponds to the approximation of 
the boundary layer of zero-width. Note the agreement with the Theorem 
2: the outer solutions agree very well. 

In Fig[8l we also compare the numerical solutions for the electric current 
in the case of the rigorous boundary conditions and in the case of the corre- 
sponding Donnan boundary condition approximation. The both solutions 
agree very well. 

6. Discussion 

The studied ID PNP model describes ion fiows in the conical geometry in 
the presence of a permanent charge on the channel wall. The dimensionality 
reduction from 3D to ID results both in an entropic potential asymmetry 
and in the appearence of an inhomogeneous, asymmetric ID charge density. 
In this paper, we presented the details of the singular perturbation approach 
to the ID PNP system including these profound effects. In another paper 
[1], these results are used further to find an analytical expression for the 
electric current and to quantify the rectification effect. When the surface 
charge density is small enough, the obtained theoretical expressions are 
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Fig. 5. (Color online) Potential profile <i>(z) for 



Don 

+ , L 



mM, eg™ L 



156 niM, p = 125 
87.3 mM, eg™ ^ = 83.6 mM, a = -0.1 e/nm^. Calculations are 

done for <P^°"{0) = -0.014 V,' ^^on^^) ^ Q_gg4 y ^j^) j^^j $Don(-Q) ^ _0.014 V, 

(jjDonj-^^-j _ —1.006 V (b). We show perturbation solution in first order (in e): 
dashed line and second order: dotted line. The numerical solution is represented 
by symbols. Insets depict closer look into left and right boundary layers. 



expected to provide good approximations for the electric potential profile, 
concentrations of ions, the total flux, as well as for the electric current. 

Although such weakly charged channels seem to be less interesting from 
the experimental point of view, they are important model systems 0]. In 
particular, based on our analytical solutions, complemented also by the 
numerical analysis, we studied and compared two different boundary condi- 
tions: the rigorous boundary conditions and the so-called Donnan bound- 
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Fig. 6. (Color online) In the upper panel (a) we show the concentration profile 
Ck+{z) for (7 = -0.02 e/nm2, and the Donnan b. c: cg^L = mM, cg^p = 105 



mM, 



„Don 

''Gl- 



ial. 8 mM, c, 



Don 
C1-. R 



103.8 mM, $Don^o) = -0.003 V, ^^on^^,) 



0.999 V; and in the lower panel (b) we show the I — U dependence. We show 
perturbation solution in first order (in e): solid line, second order: dashed line, 
third order: dotted line, and fourth order: dash-dotted line. The numerical solution 
is represented by symbols. 



ary condition. The heuristic assumption of the Donnan equilibrium at the 
boundaries of the channel [8] in a highly non-equilibrium situation (current 
rectification) should and must be questioned. However, our analysis clearly 
showed that when the channel length exceeds much the Debye screening 
length, i.e. in the limit 1/A 0, the use of the Donnan boundary condi- 
tions is well justified. 
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Fig. 7. (Color online) Potential profile $(z) for cl — cr — 0.1 M, a — —0.02 e/nm^. 
Calculations are done for $(0) = ^{L) = V (solid line) and corresponding Donnan 
b. c. $Don(Q) ^ _2.9 mV, $Don(^) ^ _2_3 jj^Y (dashed line) (a); $(0) = V, 
$(L) = IV (solid line) and corresponding Donnan b. c. $°°"(0) = -0.003 
V, $^°"(i) — 0.999 V (dashed line) (b). We show perturbation solution in first 
order (in e). Symbols in all cases stands for numerical solution, lines represent 
perturbation theory. 
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I [nA] 0.6-1 




Fig. 8. (Color online) Numerically obtained current voltage (/ — U) characteristics 
for (7 = — 0.02e/nm^. Calculation are done for: rigorous b. c. (solid line) cl = cr = 
0.1 M and corresponding Donnan b. c. (squares) c^+\ ~ ^^"^ mM, cj^+"i^ — 105 
mM, eg™ L = 131.8 niM, eg™ p = 103.8 niM. 



Appendix A 

ID model reduction 

Let us start from the 3D electro-diffusion equation 

= -V . j,{r,t) (A.l) 

where 

fiif, t) = -DiVciif, t) + ez^i/i/(r)ci(r, t) (A.2) 

defines the mass current ji of the i-th species of ions, r is the position vector, 
Di is the diffusion constant, //j is the mobility of the ion particles, and 
S{'r) = — V$(r) denotes the electric field. The consistency with the thermal 
equilibrium demands that the mobility of the ions /i, fulfills the Einstein 
relation, reading Di = fii/f3, with the inverse temperature /3 = 1//cbT. 

The electric potential $(r) is governed in a self-consistent manner by 
the Poisson equation 

eoV • [e(f)V«>(f)] = - J2 Piir) - PaA^, (A.3) 

i=K+,Cl~ 

where Pfix(^) = S{r — R{z))a, with 5 being the Dirac delta-function, repre- 
sents the fixed charges located on the inside of the channel wall. Pi{r) de- 
notes the density of mobile ions, and e(r) = epG(r — i?(z)) + e^0(— r + i?(z)) 
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describes dielectric properties of the system, with ep,e^ being the relative 
dielectric constants of the polymer and water, respectively, and being 
the Heaviside step function. In above-presented definitions we made use of 
cylindrical coordinates that are natural for the geometry given in Fig. [TJ 

The set of relations (jA.ip and ()A.3P for ji = const constitutes the system 
of coupled 3D Poisson-Nernst-Planck (PNP) equations. 

Under assumptions of instantaneous equillibration in the transverse di- 
rection as in Eq. ^ and employing the definition of divergence given below 
(App. B) we derive the set of ID PNP equations (Eqs. ([5]) and ([6])), see 
also Ref. Hdl]. 



Appendix B 

The divergence theorem. 

The definition of the divergence of a vector u is 

V • M = lim — / u ■ h dS = lim . , , , — / u ■ h dS, (B.l) 

V-^oVJs Az^o A{z)Az Js 

where V{z) = A{z)Az denotes the volume surrounded by the closed surface 
S with the normal vector h. For the flux field j{z) = j{z)ez (no ion flux 
across the channel wall) Eq. (jB.ip can be reduced [12] to the form 

Az^o A[z)Az A{z) dz 

(B.2) 

In the case of the electric field we are intrested in the divergence: V • 
(e(r, z)E{z)^ . When £{r, z) = £,{r, z)e, + f^(r, z)er the surface S is either 

parallel or perpendicular to £. Inside the channel the electric field has only 
the z component depending on z coordinate: £{z) = £z{z)ez- As the channel 
wall is charged with prescribed charge density a, we find 



(ep^pVmer «-^^) " A(" ) ^^^^^^^^^^ e^Riz)' 



(B.3) 



where A{z) = 'kB?{z). The jump in the electric field is due to the surface 
charge density (the electrostatic boundary condition on the channel wall): 

^P'^polymer ~ ^w£ = —■ (^-4) 
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Appendix C 

The outer solutions - general expressions. 

Appendix C.l Zeroth order in e 
In zeroth order in e, the solution reads 

4sU(^)=(il+c<s)ij^(.). (0.1) 

Substituting that expression for the sum of concentrations in the second 
equation, we get 

• J^Q^ = if cl = cr, (see Appendix B) 

r(0) , 



CR 



Appendix C.2 First order in e 

The solution reads 



• -^(o) - ^ 
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And the equation for the potential has the form 



(C.6) 



C(o),eI^J I J 

Appendix C.3 Higher orders in e 
The solution reads 

t(0) 



where ^-^^^ denotes the part of solution that does not depend on coeffi- 
cients of second order in e 

(z) - f _ (z) + 2:^c(°) fz) 



And the equation for the potential for n>2 has the form 
The general solution for sum of concentrations for n > 3 reads 



- / 1^ - E^SL,.!^-)^ - i»S . (CIO) 



j{0) 
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where 

(C.12) 

with 

The sohition is itcrativcly given in orders of e. Note that, the solution, 
which is vahd in the interior of the channel (called the outer expansion), 
is not valid near the both endpoints of the interval. Thus we have to find 
a different representation for the solutions near z = and z = 1. These 
are boundary layer representation, which we consider in the next section. 
The outer expansion contain unknown constants jj^^^ , l'^^^ , C'^^^ , E^^^ , where 
z = 0, 1, 2, . . ., which are determined by the matching condition. 



Appendix D 

The outer solutions - determination of constants. 

We obtain analytical formulas for the constants. In zeroth order in e, in 
the cases when 



CL = CR = c: 



■^(0) ~ ^' '-'(0) ~ 



2 



'(1) - Y^-/-R' -(1) - -^i^-R- 
• CL / cr: 

J(J) = 27r(l + 7) (cR - cl) , Cfjj = ^ (cr(1 + 7) - cl) 

r{0) _ t(0) p(0) _ ln(27r7CL) 

'(0) - ^(0) In (cl/cr) ' ""(O) - ln(cL/cR) 



(D.l) 



(D.2) 
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For higher orders these calculations get pretty quickly tedious. Moreover, 
the formulas are very long. Thus, we do not present them explicitly. Instead, 
we propose a simple algorithm that we implemented in Maple. 
Algorithm: 

. c[;]_^(z) (Eq. del])) and cl>g](z) (Eq. 



j(?) = (i + .)(tS}(o)-t{;)(i) 

c{?} = jSV,-Tf°)(0) 



.(0) 



<!>[;](.) (Eq. m) 
= solve (c^[;](o) - $[;](!) + 0[;;(-oo) - 0!;;(oo),i[;) 

(0) 

(1) 



4°}=solve(cl>[;](O)-0;;](oo),E(°) 



^{(0),...,n-l}l^i 

- for n = 2: Eq. ([091) 

- forn > 3: Eqs. ([UT2]) and I^JTL^ 



■'S = (i+7)(iS...„_„(o)-ij;!„„„_„(i) 



+ 



t(0) 

^(0) _ -^(n) ^(0) 



^{l!l...,n-i}(o)+2yp[:i(oo)n[;;(oo). 

C(n),s(^) (Eqs- dEHj) and ([OTTl) ) and $[°}(z) (Eq. ([OTOl) ) 

= -ive (ci>[S(o) - c^Sd) + 0[:;(-oo) - ^s°5(oo), ij;; 
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. 4°,' (<>(») -<!(»)'<!) 



In order to find values of |p|j^j(c>o), n^^^{oo), 0|^j(c>o) | and <^p^^^{—oo), 

fi||||(— co), (— oo) I we make use of the series expansion of tlie Donnan 

equilibrium conditions (j40p at and zr,, respectively (see Theorem 1 in 
Sec. 4). 
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